using namespace arma; 

template<class Density,class Proposal>
class Kernel_RHMCProbit<Proposal>
{
	public:
		Kernel_RHMCProbit(Density d,Proposal *p, mat S)
		{
			//_p=malloc(sizeof(Proposal));
			_p=p;
			boost::random::uniform_01<> U;
			runif=new RandomG::Random<boost::random::uniform_01<> >(U);
			_Sig=1;
			_step=10;
			eps=0.00000000000001;//_step/(4*sqrt(min(S.diag())));
			//eig_sym(eigval, eigvec, S);
			//eps=1;//sqrt(as_scalar(max(eigval)));
			cout << "eps " << eps;
			_M=inv(S);
			_d;

						
		}
		~Kernel_RHMCProbit(){
			delete runif;
		}
		void Move(mat *X, int_n)
		{
			int n=X->n_rows;//Y and X as pointers
			int m=X->n_cols;
			if(_bool==1)
			{
				Moved.zeros();
			}
			mat Y(n,m);
			Y.fill(0);
			mat xt;
			_p->Set_n(m);
			_sum=0;
			_i=0;
			for(int i=0;i<n;i++)
			{
				Y(i,span::all)=Prop((*X)(i,span::all));//Y[i]=p->r(X[i])
				mat xt=Y.row(i);
				xt=Accept(X->row(i),Y.row(i),_n);
				(*X)(i,span::all)=xt;			
				_i++;
			}
			_accept=(double)_sum/n;

			cout << "accept = " << (double)_sum/n;

		}
		mat Prop(mat X)
		{
			mat Y=X;
			mat ps;
			mat Xs;
			mat p;
			_Y1=Y;
			double th=0.00000001;
			for(int i=0;i<_step;i++)
			{	
				mat p1=p;
				while(err>th)
				{
					ps=f1(p,X,eps,p1);
					err=norm(p,ps);
					p=ps;
				}
				mat X1=Y;
				while(err>th)
				{
					Xs=f2(p,X,eps,X1);
					err=norm(X,Xs);
					p=ps;
				}
				p=f3(p,X,eps);

			}
			Y=X;
			return Y;

		}
		mat f1(mat p,mat X,double eps,mat p1)
		{
			return p1-0.5*eps*Ht(X,p);			
		}
		mat f2(mat p,mat X,double eps,mat X1)
		{
			return X1+0.5*eps*(Hp(X,p)+Hp(X1,p));			
		}
		mat f3(mat p, mat X, double eps)
		{
			return p-0.5*eps*Ht(X,p);
		}
		mat Ht(mat X, mat p)
		{

		}
		mat Hp(mat X, mat p)
		{

		}
		mat InfoFisher(mat X)
		{
		}
		mat Grad_info(mat X,int i)
		{
		}



	private:
		double eps;
		double _Sig;
		int _Step;
		RandomG::Random<boost::random::uniform_01<> > *runif;
		Proposal _p;
		Density _d;
		mat _M;

		mat _Y1;

};
